Biostatistiques

Serge-Étienne Parent

09 février 2019

Objectifs spécifiques (1/2)

À la fin de ce chapitre, vous

  • serez en mesure de définir les concpets de base en statistique: population, échantillon, variable, probabilité et distribution
  • serez en mesure de calculer des statistiques descriptives de base: moyenne et écart-type, quartiles, maximum et minimum
  • comprendrez les notions de test d’hypothèse, d’effet et de p-value, ainsi qu’éviter les erreurs communes dans leur interprétation

Objectifs spécifiques (2/2)

  • saurez effectuer une modélisation statistique linéaire simple, multiple et mixte, entre autre sur des catégories
  • saurez effectuer une modélisation statistique non linéaire simple, multiple et mixte

Tout commence par une question

Que voulez-vous exprimer?

Source

La question

Que voulez-vous exprimer?

Source

Où renforcer l’avion?

Disposition des trous de balle sur les avions britaniques de retour de leur sortie durant la deuxième guerre mondiale.

Source: Whitlock et Schuluter (2015), The Analysis of Biological Data.

Définitions

La statistique est l’étude des méthodes pour mesurer des aspects de populations à partir d’échantillons et pour quantifier l’incertitude des mesures. (ma traduction de Whitlock et Schuluter (2015)).

Population

Ensemble circonscrit du sujet d’étude.

Échantillon

Sous-ensemble de la population.

Si l’échantillon est représentatif de la population (non biaisé), on peut extrapoler ses propriétés à la population (inférence).

Variable aléatoire

Dimension d’un objet catractérisé par une distribution de probabilité.

Probabilité

Une probabilité est la vraisemblance qu’un évènements se réalise chez un échantillon.

Distribution de probabilité

Une distribution décrit la probabilité d’obtenir une valeur (distribution discrète) ou une plage de valeurs (distribution continue) dans une échantillon pris au hasard d’une population.

Une distribution des probabilité est décrite par un ou plusieurs paramètres, par exemple une moyenne et un écart-type dans le cas d’une distribution normale.

Une statistique

Une statistique est une estimation d’un paramètre calculée à partir des données, par exemple une moyenne et un écart-type échantillonnaux.

Statistiques fréquentielles ou bayésiennes

Fréquentielle. Les données sont générées par des mécanismes stochastiques décrits par des distributions de probabilités, et nous cherchons à déterminer la probabilité que les données soient générées par ces mécanismes (p-valule). <– approche la plus commune (et couverte dan ce chapitre)

Statistiques fréquentielles ou bayésiennes

Bayésienne. Les paramètres et leur incertitude sont déterminés par les données et les connaissances préalables. <– approche de plus en plus utilisée (effleurée dans le chapitre 6, en extra)

Bayes ou freq?

Tout dépend de la question posée.

  • Bayésien: quelle est la probabilité qu’il y ait de la vie sur Mars?
  • Fréquentiel: est-ce que les données sont conformes avec avec l’hypothèse de la vie sur Mars

exemple tirée du blogue Dynamic Ecology

Plusieurs types de distribution

Les tests d’hypothèse

Un test d’hypothèse permet de décider si une hypothèse est confirmée ou rejetée à un seuil de probabilité prédéterminé.

H0. L’hypothèse nulle propose l’absence d’effet statistique (c’est l’hypothèse de l’avocat du diable 😈) .

H1. L’hypothèse alternative propose l’inverse (😇).

Exemple

Est-ce que les données sont conforment avec avec l’hypothèse de la vie sur Mars.

H0. Il n’y a pas de vie sur Mars.

H1. Il y a de la vie sur Mars.

Dans une expérience, la confirmation de l’hypothèse nulle n’est pas un échec.

Test de t (Student) à deux échantillons

Comparaison des moyennes deux échantillons dont la distribution est normale et dont la variance est la même.

## # A tibble: 5 x 2
##   echantillon value
##   <chr>       <dbl>
## 1 a           10.1 
## 2 b            7.98
## 3 a            5.82
## 4 b           11.8 
## 5 a            9.55
## 
##  Two Sample t-test
## 
## data:  value by echantillon
## t = 1.5581, df = 58, p-value = 0.1247
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2990577  2.3995879
## sample estimates:
## mean in group a mean in group b 
##        9.296136        8.245871

Enregistrer les résultats dans un objet

## List of 9
##  $ statistic  : Named num 1.56
##   ..- attr(*, "names")= chr "t"
##  $ parameter  : Named num 58
##   ..- attr(*, "names")= chr "df"
##  $ p.value    : num 0.125
##  $ conf.int   : num [1:2] -0.299 2.4
##   ..- attr(*, "conf.level")= num 0.95
##  $ estimate   : Named num [1:2] 9.3 8.25
##   ..- attr(*, "names")= chr [1:2] "mean in group a" "mean in group b"
##  $ null.value : Named num 0
##   ..- attr(*, "names")= chr "difference in means"
##  $ alternative: chr "two.sided"
##  $ method     : chr " Two Sample t-test"
##  $ data.name  : chr "value by echantillon"
##  - attr(*, "class")= chr "htest"
## [1] 0.1246579

p-value (1/)

“La p-value n’a jamais été conçue comme substitut au raisonnement scientifique” Ron Wasserstein, directeur de l’American Statistical Association [ma traduction].

p-value (2/)

Un résultat montrant une p-value plus élevée que 0.05 est-il pertinent?

Lors d’une conférence, Dr Evil ne présentent que les résultats significatifs de ses essais au seuil de 0.05. Certains essais ne sont pas significatifs, mais bon, ceux-ci ne sont pas importants…

Rappel

Dans une expérience, la confirmation de l’hypothèse nulle n’est pas un échec.

p-value: erreurs (3/)

  1. Pour évaluer l’importance de l’effet: voir le coefficient. Pour évaluer son incertitude: voir son intervalle de confiance.
  2. Il est tout aussi important de savoir que le traitement fonctionne que de savoir qu’il ne fonctionne pas.
  3. Le seuil de 0.05 est arbitraire.

p-value: p-hacking ou data dregging (4/)

Modifier les données ou la question posée dans le but d’obtenir une p-value < 0.05, c’est tricher.

L’analyse de variance (ANOVA)

Comparaison des moyennes de plus de 2 groupes (à variance égale)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## Species       2  437.1  218.55    1180 <2e-16 ***
## Residuals   147   27.2    0.19                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modélisation statistique

\[Y = f\left( X \right)\] \(X\): variables prédictives / indépendantes / covariables / explicatives \(Y\): variables de sortie / réponse / dépendante

Modèles prédictifs / explicatifs

⚠️ Ne pas confondre. ⚠

Modèles prédictifs. Le modèle doit être testé sur des données qui n’ont pas servi à l’entraîner. (chapitre 12)

Modèle explicatif. Le modèle sert à tester des hypothèses. (chapitres 5)

Variables fixes et aléatoires

Variables fixes: testées lors de l’expérience: dose du traitement, espèce/cultivar, météo, etc.

Variables aléatoires: sources de variation qui génèrent du bruit dans le modèle: les unités expérimentales ou le temps lors de mesures répétées.

Modèle linéaire univarié

\[y = \beta_0 + \beta_1 x + \epsilon\]

Modèle linéaire univarié

Question. Le rendement varie-t-il selon la dose d’azote?

H0. la dose d’azote n’affecte pas le rendement (\(\beta_1 = 0\)).

La fonction lm

## 
## Call:
## lm(formula = yield ~ nitro, data = lasrosas.corn)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.183 -15.341  -3.079  13.725  45.897 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 65.843213   0.608573 108.193  < 2e-16 ***
## nitro        0.061717   0.007868   7.845 5.75e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.66 on 3441 degrees of freedom
## Multiple R-squared:  0.01757,    Adjusted R-squared:  0.01728 
## F-statistic: 61.54 on 1 and 3441 DF,  p-value: 5.754e-15

predict

Analyse des résidus (1/)

Résidus: erreurs du modèle, \(\epsilon = y - \hat{y}\), fonction residuals().

##      nitro residus_lm
## 3086  99.8  -14.27259
## 231  106.0  -16.50523
## 2911 124.6  -35.33317
## 800  106.0  -18.80523

Analyse des résidus (2/)

  • Pas de structure identifiable dans les résidus
  • Distribution normale

Modélisation linéaire multiple

\[ y = X \beta + \epsilon \]

Modèles linéaires univariés avec variable catégorielle nominale

topo est la classe topographique

##    year       lat      long yield nitro topo     bv rep nf
## 1  1999 -33.05113 -63.84886 72.14 131.5    W 162.60  R1 N5
## 35 1999 -33.05155 -63.84645 56.70 131.5   HT 169.49  R1 N5
## 51 1999 -33.05174 -63.84532 59.94 131.5    E 182.12  R1 N5
## 67 1999 -33.05195 -63.84412 70.35 131.5   LO 168.00  R1 N5

Matrice-modèle avec intercept

##    (Intercept) topoHT topoLO topoW
## 1            1      0      0     1
## 35           1      1      0     0
## 51           1      0      0     0
## 67           1      0      1     0

Matrice-modèle sans intercept

##    topoE topoHT topoLO topoW
## 1      0      0      0     1
## 35     0      1      0     0
## 51     1      0      0     0
## 67     0      0      1     0

Modèle linéaire sur des catégories = régression multiple

## 
## Call:
## lm(formula = yield ~ topo, data = lasrosas.corn)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.371 -11.933  -1.593  11.080  44.119 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  78.6653     0.5399 145.707   <2e-16 ***
## topoHT      -30.0526     0.7500 -40.069   <2e-16 ***
## topoLO        6.2832     0.7293   8.615   <2e-16 ***
## topoW       -11.8841     0.7039 -16.883   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.59 on 3439 degrees of freedom
## Multiple R-squared:  0.4596, Adjusted R-squared:  0.4591 
## F-statistic:   975 on 3 and 3439 DF,  p-value: < 2.2e-16

Du modèle linéaire à l’anova

## Analysis of Variance Table
## 
## Response: yield
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## topo         3 622351  207450  974.96 < 2.2e-16 ***
## Residuals 3439 731746     213                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modèles linéaires univariés avec variable catégorielle ordinale

##   (Intercept) statut_o.L statut_o.Q    statut_o.C statut_o^4
## 1           1 -0.6324555  0.5345225 -3.162278e-01  0.1195229
## 2           1 -0.3162278 -0.2672612  6.324555e-01 -0.4780914
## 3           1  0.0000000 -0.5345225 -4.095972e-16  0.7171372
## 4           1  0.3162278 -0.2672612 -6.324555e-01 -0.4780914
## 5           1  0.6324555  0.5345225  3.162278e-01  0.1195229
## attr(,"assign")
## [1] 0 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$statut_o
## [1] "contr.poly"

Régression linéaire multiple

## 
## Call:
## lm(formula = yield ~ lat + long + nitro + topo + bv, data = lasrosas.corn)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.405 -11.071  -1.251  10.592  40.078 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.946e+05  3.309e+04   5.882 4.45e-09 ***
## lat          5.541e+03  4.555e+02  12.163  < 2e-16 ***
## long         1.776e+02  4.491e+02   0.395    0.693    
## nitro        6.867e-02  5.451e-03  12.597  < 2e-16 ***
## topoHT      -2.665e+01  1.087e+00 -24.520  < 2e-16 ***
## topoLO       5.565e+00  1.035e+00   5.378 8.03e-08 ***
## topoW       -1.465e+01  1.655e+00  -8.849  < 2e-16 ***
## bv          -5.089e-01  3.069e-02 -16.578  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.47 on 3435 degrees of freedom
## Multiple R-squared:  0.5397, Adjusted R-squared:  0.5387 
## F-statistic: 575.3 on 7 and 3435 DF,  p-value: < 2.2e-16

Standardiser les données pour comparer les effets

## 
## Call:
## lm(formula = yield ~ lat + long + nitro + topo + bv, data = lasrosas.corn_sc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.405 -11.071  -1.251  10.592  40.078 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  78.9114     0.6666 118.376  < 2e-16 ***
## lat           3.9201     0.3223  12.163  < 2e-16 ***
## long          0.3479     0.8796   0.395    0.693    
## nitro         2.9252     0.2322  12.597  < 2e-16 ***
## topoHT      -26.6487     1.0868 -24.520  < 2e-16 ***
## topoLO        5.5647     1.0347   5.378 8.03e-08 ***
## topoW       -14.6487     1.6555  -8.849  < 2e-16 ***
## bv           -4.9253     0.2971 -16.578  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.47 on 3435 degrees of freedom
## Multiple R-squared:  0.5397, Adjusted R-squared:  0.5387 
## F-statistic: 575.3 on 7 and 3435 DF,  p-value: < 2.2e-16

Intéraction

Une intéraction est une pente additionnelle informant sur l’effet statistique combiné de deux variables.

Interface-formule: le symbole : appelle l’intéraction et le synbole * appelle les effets simples et les intéractions.

## 
## Call:
## lm(formula = yield ~ nitro * topo, data = lasrosas.corn_sc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.984 -11.985  -1.388  10.339  40.636 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   78.6999     0.5322 147.870  < 2e-16 ***
## nitro          1.8131     0.5351   3.388 0.000711 ***
## topoHT       -30.0052     0.7394 -40.578  < 2e-16 ***
## topoLO         6.2026     0.7190   8.627  < 2e-16 ***
## topoW        -11.9628     0.6939 -17.240  < 2e-16 ***
## nitro:topoHT   1.2553     0.7461   1.682 0.092565 .  
## nitro:topoLO   0.5695     0.7186   0.792 0.428141    
## nitro:topoW    0.7702     0.6944   1.109 0.267460    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.38 on 3435 degrees of freedom
## Multiple R-squared:  0.4756, Adjusted R-squared:  0.4746 
## F-statistic: 445.1 on 7 and 3435 DF,  p-value: < 2.2e-16

Modèles linéaires généralisées

Évaluer l’effet sur une variable entière, booléenne, factorielle, etc.

La fonction glm

## 
## Call:
## glm(formula = worms ~ trt, family = "poisson", data = cochran.wireworms)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8279  -0.9455  -0.2862   0.6916   3.1888  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.1823     0.4082   0.447 0.655160    
## trtM          1.6422     0.4460   3.682 0.000231 ***
## trtN          1.7636     0.4418   3.991 6.57e-05 ***
## trtO          1.5755     0.4485   3.513 0.000443 ***
## trtP          1.3437     0.4584   2.931 0.003375 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 64.555  on 24  degrees of freedom
## Residual deviance: 38.026  on 20  degrees of freedom
## AIC: 125.64
## 
## Number of Fisher Scoring iterations: 5

Sur graphique

Modèles non-linéaires

Mitscherlich

\[ y = A \left( 1 - e^{-R \left( E + x \right)} \right) \]

La focntion nls

## 
## Formula: yield ~ A * (1 - exp(-R * (E + nitro)))
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## A 75.023427   3.331860  22.517   <2e-16 ***
## E 66.164111  27.251592   2.428   0.0184 *  
## R  0.012565   0.004881   2.574   0.0127 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.34 on 57 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 8.063e-06

Résidus

Modèles à effets mixtes

À la différence d’un effet fixe, un effet aléatoire (\(b\)) sur des variables aléatoires (\(Z\)) sera toujours distribué normalement avec une moyenne de 0 et une certaine variance. Pour un modèle à intercept aléatoire, nous aurons

\[ y = X \beta + Z b + \epsilon \]

Pour un modèle à pente aléatoire, nous aurons (à vérifier)

\[ y = X \left( \beta + Zb \right) + \epsilon \]

La fonction lme du module nlme

random = ~ slope|intercept

## Linear mixed-effects model fit by REML
##  Data: lasrosas.corn 
##        AIC      BIC    logLik
##   26535.37 26602.93 -13256.69
## 
## Random effects:
##  Formula: ~1 | year
##         (Intercept)
## StdDev:    20.35425
## 
##  Formula: ~1 | rep %in% year
##         (Intercept) Residual
## StdDev:    11.17447 11.35617
## 
## Fixed effects: yield ~ lat + long + nitro + topo + bv 
##                  Value Std.Error   DF    t-value p-value
## (Intercept) -1379436.9  55894.55 3430 -24.679273   0.000
## lat           -25453.0   1016.53 3430 -25.039084   0.000
## long           -8432.3    466.05 3430 -18.092988   0.000
## nitro              0.0      0.00 3430   1.739757   0.082
## topoHT           -27.7      0.92 3430 -30.122438   0.000
## topoLO             6.8      0.88 3430   7.804733   0.000
## topoW            -16.7      1.40 3430 -11.944793   0.000
## bv                -0.5      0.03 3430 -19.242424   0.000
##  Correlation: 
##        (Intr) lat    long   nitro  topoHT topoLO topoW 
## lat     0.897                                          
## long    0.866  0.555                                   
## nitro   0.366  0.391  0.247                            
## topoHT  0.300 -0.017  0.582  0.024                     
## topoLO -0.334 -0.006 -0.621 -0.038 -0.358              
## topoW   0.403 -0.004  0.762  0.027  0.802 -0.545       
## bv     -0.121 -0.012 -0.214 -0.023 -0.467  0.346 -0.266
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.32360269 -0.66781575 -0.07450856  0.61587533  3.96434001 
## 
## Number of Observations: 3443
## Number of Groups: 
##          year rep %in% year 
##             2             6

Effet aléatoire

## Level: year 
##      (Intercept)
## 1999   -13.71488
## 2001    13.71488
## 
## Level: rep %in% year 
##         (Intercept)
## 1999/R1 -12.3828119
## 1999/R2  -2.0608423
## 1999/R3  10.3099880
## 2001/R1  -9.2956080
## 2001/R2   0.8470003
## 2001/R3  12.5822738
##   mean(`(Intercept)`)
## 1       -2.279475e-11

Modèles mixtes non-linéaires

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: yield ~ A * (1 - exp(-R * (E + nitro))) 
##  Data: engelstad.nitro 
##        AIC      BIC    logLik
##   477.2285 491.8889 -231.6143
## 
## Random effects:
##  Formula: A ~ year | loc
##  Structure: General positive-definite, Log-Cholesky parametrization
##               StdDev       Corr  
## A.(Intercept)  2.602881636 A.(In)
## A.year         0.003066058 -0.556
## Residual      11.152760033       
## 
## Fixed effects: list(A ~ 1, E ~ 1, R ~ 1) 
##                  Value Std.Error DF   t-value p-value
## A.(Intercept) 74.58277  4.722806 56 15.792046  0.0000
## E             65.57006 25.534747 56  2.567876  0.0129
## R              0.01308  0.004807 56  2.720245  0.0087
##  Correlation: 
##   A.(In) E     
## E  0.379       
## R -0.483 -0.934
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.83376801 -0.89290179  0.07419027  0.68354461  1.82431474 
## 
## Number of Observations: 60
## Number of Groups: 2

Modèle mixte non-linéaire

Modélisation non-linéaire multiple

–> Parent et al., 2018

Objectifs spécifiques (1/2)

À la fin de ce chapitre, vous

  • serez en mesure de définir les concpets de base en statistique: population, échantillon, variable, probabilité et distribution
  • serez en mesure de calculer des statistiques descriptives de base: moyenne et écart-type, quartiles, maximum et minimum
  • comprendrez les notions de test d’hypothèse, d’effet et de p-value, ainsi qu’éviter les erreurs communes dans leur interprétation

Objectifs spécifiques (2/2)

  • saurez effectuer une modélisation statistique linéaire simple, multiple et mixte, entre autre sur des catégories
  • saurez effectuer une modélisation statistique non linéaire simple, multiple et mixte